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Abstract 

A traction-displacement relationship that may be embedded into a cohesive zone model for 
microscale problems of intergranular fracture is extracted from atomistic molecular-dynamics 
simulations. A molecular-dynamics model for crack propagation under steady-state conditions is 
developed to analyze intergranular fracture along a flat 299 [110] symmetric tilt grain boundary 
in aluminum. Under hydrostatic tensile load, the simulation reveals asymmetric crack 
propagation in the two opposite directions along the grain boundary. In one direction, the crack 
propagates in a brittle manner by cleavage with very little or no dislocation emission, and in the 
other direction, the propagation is ductile through the mechanism of deformation twinning. This 
behavior is consistent with the Rice criterion for cleavage vs. dislocation blunting transition at 
the crack tip. The preference for twinning to dislocation slip is in agreement with the predictions 
of the Tadmor and Hai criterion. A comparison with finite element calculations shows that while 
the stress field around the brittle crack tip follows the expected elastic solution for the given 
boundary conditions of the model, the stress field around the twinning crack tip has a strong 
plastic contribution. Through the definition of a Cohesive-Zone- Volume-Element - an atomistic 
analog to a continuum cohesive zone model element - the results from the molecular-dynamics 
simulation are recast to obtain an average continuum traction-displacement relationship to 
represent cohesive zone interaction along a characteristic length of the grain boundary interface 
for the cases of ductile and brittle decohesion. 
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1. Introduction 

Cohesive zone models (CZMs) approximate traction-displacement relationships along an 
interface (Tvergaard and Hutchinson, 1992; Costanzo and Allen, 1995; Camacho and Ortiz, 
1996; Klein and Gao, 1998) and are frequently used in conjunction with the finite element 
method (FEM) to study fracture in a wide variety of materials. Idealized traction-displacement 
behavior of interface debonding is embedded into CZM elements. In fracture studies at the 
micromechanical level, CZM elements can be placed between the continuum finite elements that 
discretize the grain interior to predict transgranular fracture or placed between the continuum 
finite elements on either side of a grain boundary to predict intergranular fracture. 

Modeling of material failure with CZMs has been advanced to the level of being able to 
perform large scale simulations of fracture in polycrystals. Recently, Zavattieri et al. (2001) 
studied the fracture of alumina-ceramic microstructures subjected to multi-axial dynamic 
loading. The effective size of the polycrystalline specimen studied was 0.54 x 0.19 mm, thus 
reaching macroscopic scales. A bilinear traction-displacement relationship parameterized to 
empirical data, such as macroscopic fracture toughness K f c was used. Zavattieri and Espinosa 
(2003) used a modification of the same model to study interface effects of an alumina specimen 
in contact with steel plates. Wei and Anand (2004) have used a modified CZM model to study 
intergranular fracture in nanocrystalline Ni. In their finite element (FE) model simulation, the 
CZM element approximated both reversible and irreversible inelastic sliding-separation 
deformations at the grain boundaries prior to failure. The parameterization of the model was, 
again, performed by using available experimental data for stress-strain curves of nanocrystalline 
Ni in tension with an average grain size of 15-40 mn and having a comparatively large number 
of grains. Iesulauro et al. (2002) have applied the CZM technique to simulate fatigue crack 
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initiation in A1 polycrystals. The use of MD simulations to parameterize the traction- 
displacement curve was suggested, but the actual parameterization was performed by using the 
established macroscale yield properties of aluminum. 

The macroscale values of strength and toughness that are input to the CZM in these 
references represent the aggregate responses of thousands or millions of grains, grain boundaries, 
and defects within the specimens from which they were obtained. Thus, these macroscale values 
do not represent the unique response of a particular interface at which a local fracture event 
might occur. If the microscale predictions are to become quantitative, consideration of the local 
nanoscale properties is required. One possible means of making this connection is to use the 
results of atomistic MD models as input to the CZM. This connection would allow more realistic 
simulations leading to accurate predictions of the failure properties of a large class of materials 
and microstructures, even when experimental data is not available. 

Attempts to extract relevant parameters for the decohesion law of a CZM from atomistic 
(molecular-dynamics or molecular-static) simulations have been made by various groups in the 
last few years (Gall et ah, 2000; Komanduri et al., 2001; Spearot et ah, 2004). The approach in 
all of these works is based on simulating the debonding of a flat interface under a constant tensile 
strain rate perpendicular to the interface. In these references, the system size is between 4 and 8 
mn, and the dynamics of the atoms is severely constrained by the boundary conditions, which do 
not allow for Poisson lateral contraction and shear deformation. As a result, plastic processes, 
such as dislocation slip, are strongly suppressed. Consequently, the simulated mechanism for 
interface decohesion in these references reproduces the process of atomic adhesion (strength) 
rather than that of fracture at the interface. Raynolds et al. (1996) used a similar setup to study 
adhesion in an NiAl-Cr interface by first principles calculations. 
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The boundary conditions, at which a typical CZM element operates in a large scale FE model 
are very different from the ones used in the referenced MD and first principles simulations. 
Typically, the CZM elements are embedded within a system of finite elements, that reproduce 
the elastic and plastic response of the surrounding material to both the external load and the 
crack-tip stress. In contrast, decohesion parameters, such as peak stress and opening 
displacement of the CZM curve, are extracted from an atomistic volume less than 10 mn in each 
direction. The lack of an adequate surrounding volume of material suppresses the plastic 
processes, such as dislocation nucleation, limiting the accommodation of defonnation at the 
interface and forcing it to debond in an unnatural manner. The periodic boundary conditions 
usually applied in these models cause the simulation to create a response of an array of repeating 
units with a strong overlap of image elastic forces rather than the response of a single specimen 
unit. Consequently, the resulting decohesion curves cannot be directly applied to derive the 
constitutive laws for CZM elements. 

The main goal of the approach described in the present study is to extract, and understand the 
contributions to, an MD-based CZM decohesion law for intergranular fracture under local 
conditions similar to those experienced by the CZM element in a polycrystalline FE model. The 
CZM decohesion law reflects the response of the CZM element to an approaching and 
propagating crack (Costanzo and Allen, 1995; Davila, 2001). Thus, the MD model should be a 
model of crack propagation rather than of adhesion. The MD model used in this study is built to 
simulate a crack propagating through a flat high-energy grain-boundary in aluminum (Yamakov 
et ah, 2005). 

The paper is constructed as follows: The simulation approach is described in Section 2. The 
mechanism of intergranular crack growth together with the plastic processes near the crack tip as 
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revealed by MD simulations at the atomic level is broadly discussed in Section 3. The 
contribution of the plastic processes to the stress field near the crack tip is discussed in Section 4. 
Section 5 describes a methodology for extracting a constitutive relation for a continuum CZM 
element from the MD results. The main conclusions of this study are outlined in Section 6. 


2. The simulation approach 

The simulation approach used in this study is based on a MD simulation model of crack 
propagation under time-independent, or steady-state, conditions through a flat grain boundary 
(GB) in A1 (modeled by the interatomic potential of Mishin et al. (1999)) at low temperature 
(100 K). The purpose of the MD simulation is to reveal and analyze the atomistic processes 
taking place near the crack tip and to derive a statistical traction-displacement relationship for a 
continuum CZM element. The simulation must also provide a study on the influence of the 
atomistic processes on the resulting traction-displacement (decohesion) curve. 

To produce reliable and time-independent statistics for extracting the CZM decohesion law 
the simulation model is based on a well known continuum-elastic analytical model for steady- 
state crack propagation, which is briefly reviewed in Section 2.1. Following the analytical model, 
the MD set up and the details of the MD simulation are discussed in Section 2.2. Major attention 
is given to the crystallography of the atomistic MD system, which determines the behavior of the 
model and, to a large extent, predefines the ongoing atomistic processes. 

The values of the initial and the simulation parameters of the MD system are crucial for 
recovering the steady-state regime in the MD simulation. To set up these parameters, specifically 
the system size, the initial crack length and the initial loading conditions, a continuum elastic 
solution of the system is needed. This solution cannot be taken from the analytical model, which 
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has an infinite dimension (see Section 2.1), because the finite size effects always present in an 
MD simulation make it impossible for the MD system to reproduce exactly the analytical model. 
In this work, a linear elastic solution is obtained by performing a FE simulation on an elastically 
equivalent continuum system (described in Section 2.3), and provides for a direct comparison 
between the MD and the FE results. Such a comparison helps to distinguish the role of the plastic 
processes such as dislocation nucleation, vacancy and void formation, etc. present in the MD 
simulation, but not in the FE simulation. 

2.1 The analytical model 

The theoretical model for steady-state crack propagation adapted and used in this work is 
discussed and analyzed by Langer and co-authors (Barber et al., 1989; Langer, 1992; Langer, 
1993; Langer and Nakanishi, 1993). The model represents a laterally strained strip of elastic 
material. The strip is infinite in the longitudinal direction, and a semi-infinite crack propagates 
through it (Barber et al., 1989). Far ahead of the crack tip, the strip is under uniform tension. Far 
behind the tip, the tension is relieved, and the crack opening becomes constant. Thus, there are 
two stable states for the material: “closed”, in front of the crack tip, and “open”, behind the crack 
tip. The crack propagation can be viewed as a continuous steady-state transition between these 
two states (Barber et al., 1989). The stress field of a steadily moving crack for this model was 
calculated by Ching (1994). Applying this model as a base for the MD simulation (Yamakov et 
al., 2005) ensures that after some initial unsteady growth, the crack propagation will proceed 
under steady-state conditions independent of the crack length. 

2.2 The molecular-dynamics model 
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The MD simulation is based on the previously described analytical model for steady-state 
crack propagation and has the configuration shown in Fig. 1 (Yamakov et ah, 2005). With 
periodic boundary conditions in all directions, the model represents an aluminum multilayer 
system of alternating sets of thick and thin crystalline layers separated by four flat GBs. The two 
broad layers, marked as “Crystal I” and “Crystal II”, form a bicrystalline system with a flat GB 
in the middle, through which the crack propagates. The crystallographic orientations of Crystal I 
and Crystal II are presented in Fig. 2. In the imposed coordinate system of the model, the 

orientation of Crystal I is: (x:[7 7 10]; y:[5 5 7]; z: [1 1 0]), and the orientation of 

Crystal II is: (x:[7 7 10]; y:[5 5 7]; z: [1 1 0]) (see Fig. 2). In this way, Crystal II is a 
mirror image of Crystal I relative to the crystallographic plane {5 5 7 } , which becomes the 

plane of the GB between them. The GB thus fonned is a <1 1 0> £99 symmetric tilt GB, for 
which the atomic structure in A1 is known from the literature (Dahmen, 1990). This is a high- 
angle grain boundary (tilt angle of 89.42°) with a large excess (i.e., above the perfect crystal) 
energy, y G B = 0.60±0.05 J/m 2 , estimated here for a relaxed structure at T = 100 K. The high 
excess GB energy facilitates its decohesion (Wolf, 1990). The surface energy of the GB plane y s 
at 100 K is estimated in this work at y s = 0.952±0.010 J/m 2 and is in good agreement with 
experimental data for the A1 surface energy. 

The two smaller layers, Absorbing Layer I and Absorbing Layer II, on both sides of the 
bicrystalline system (Fig. 1) have the same crystallographic orientations as Crystal II and Crystal 
I, respectively. Consequently, the GBs formed by these layers are of the same crystallographic 
type as the GB between Crystal I and Crystal II. The purpose of these layers in the simulation is 
to serve as shock- wave absorbers (Gumbsch, 1997), where a damping friction coefficient is 
applied to the atoms to absorb the phonon waves generated from the crack tips. In addition, the 
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GBs between these layers and Crystal I and II act as absorbers for the dislocations that may be 
emitted and spread out from the crack tip during propagation. In this way, the negative effect of 
the periodic boundary conditions in the y-direction creating periodic images of all crack-tip 
disturbances and influencing the crack propagation is suppressed. 

The interaction between individual atoms in the system is presented by a many-body 
embedded-atom method (EAM) potential of Mishin et al. (1999) fitted to give the correct zero- 
temperature lattice constant, a 0 = 4.05 A, elastic constants, cohesive energy, vacancy formation 
energy, etc. Of particular importance for the simulation of fracture and dislocation plasticity is 
the close fit of the potential to the experimentally measured surface and stacking-fault energies. 
Potential-dependent parameters that will be needed in this study are the relaxed stable stacking- 
fault energy, y s f = 0.146 J/m 2 , the unstable stacking-fault energy, y us = 0.168 J/m 2 , as defined in 
Mishin et al. (1999), and the unstable twinning energy, y ut = 0.210±0.010 J/m 2 estimated here 
according to the method described in Tadmor and Hai (2003). 

The system thickness h in the z-direction equals only 10(2 2 0) crystallographic planes 
(accounting for the symmetry of the fee lattice), or h = 10V2/2a o ~ 2.9 mn. This thickness is more 
than four times larger than the range of the interatomic potential, r c = 1.55a 0 = 0.628 mn (Mishin 
et al., 1999), which prevents interference of the atoms with their periodic images and preserves 
the local three-dimensional (3D) physics in the system. The small thickness in the z-direction 
allows the system size in the x- and y- directions to extend up to 21 <7 7 10>a o ~ 120 mn, and 
24<5 5 7>a 0 ~ 97 mn, respectively (see Fig 1 and Fig. 2), while limiting the number of simulated 
atoms to 1,994,000, allowing the simulation to be carried out on a modest Beowulf cluster. 

The choice of the [110] crystallographic orientation of the smallest system dimension z is 

not random (Yamakov et al. 2001). The [1 1 0] direction is the common axis of the (1 1 1) and 
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(1 1 1) planes, which are glide planes for the most common slip dislocations (i.e., 1/2<110> 
dislocations) in fee metals. It has been shown (Yamakov et ah, 2003a) that all of the six slip 
systems available on these two planes can operate unobstructed by the small system size in this 
direction. Moreover, a number of complex dislocation interactions and events are still possible 
between these six slip systems, as in a full 3D space (Yamakov et ah, 2002; 2003a). As the 
dislocation activity is expected to be very important in this study, the choice of the [110] 
direction as a columnar axis for this quasi-two dimensional set-up ensures the best resemblance 
to a full 3D environment, the main constraint being that the dislocation lines of all possible 
dislocations have to be straight lines parallel to the columnar direction. A comparison between 
this quasi-2D and a full 3D environment in a nanocrystalline model shows that some grain size 
effects, related to the curvature of the dislocation loops presented in 3D, are suppressed 
(Yamakov et ah, 2003b). Specifically, as the stress needed to expand a dislocation loop is 
inversely proportional to the loop radius, which, by itself, is limited by the very small grain size 
of the nanocrystalline metal, there is a substantial grain size dependence of the yield stress. This 
geometrical dependence is absent in a quasi-2D model, because the small thickness of the system 
constrains the dislocation lines to be always straight and parallel to each other preventing the 
formation of dislocation loops. The lack of curvature of the dislocation lines, as well as of the 
crack front, should also be taken into account as a simplification in this model. In the case where 
the grain size is large enough not to be a governing parameter, the presence of curvature usually 
helps the defect nucleation process at or around the crack tip, such as dislocations and 
microvoids. While it is expected that neglecting the effects of curvature would not qualitatively 
change the fracture mechanism, it may affect the process quantitatively in terms of slightly 
decreased peak stress and work of debonding. 
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The simulations are run at a temperature of T = 1 00 K, which is low enough to suppress GB 
and surface diffusion processes and to facilitate brittle fracture in aluminum. The temperature is 
kept constant by using the Nose-Hoover thermostat (Nose, 1984). 

The system is driven to the initial condition of the simulation before crack nucleation in 
several steps. First, the system is constructed out of four blocks of perfect crystals (each one for 
Crystal I, Crystal II, Absorbing Layer I, and Absorbing Layer II in Fig. 1), which have been 
independently thermally equilibrated at 100 K at zero constant pressure. The constant pressure 
equilibration is achieved by using the Parrinello-Rahman constant-stress simulation (Parrinello 
and Rahman, 1981). When the perfect crystals are joined together into a multilayer system to 
form 299 GBs, they are also allowed to shift along the GB planes (in the x- direction) to obtain 
the lowest energy initial GB configuration. After being assembled, the whole system is 
equilibrated at 100 K and zero pressure to achieve thermal equilibrium of the GB interfaces. 
Because the building crystal blocks have been already preheated to 100 K, the thermalization of 
the whole system takes much less time than if it were built out of ideal crystals at 0 K and then 
thermalized. 

After the thermal equilibration at zero pressure, the system is loaded hydrostatically in 
tension, i.e., 

o xx =o yy =o zz =o ( 1 ) 

and is dynamically equilibrated at this constant stress. After establishing equilibrium between the 
strain in the system and the applied external stress, the system size in all three dimensions is 
fixed under the constant strained condition. 

The transition from a constant stress to a constant strain simulation transforms the volume 
fluctuations, always present in a finite system under thermodynamic equilibrium, into stress 


10 



To be submitted to: Journal of the Mechanics and Physics of Solids 


fluctuations. Thus, further equilibration at constant strain is necessary to smooth out these 
fluctuations. The simulation then proceeds under this constant strain - constant temperature 
condition. Although the simulations are carried out under constant strain, for convenience in 
presentation, the various analyses will reference the value of prestresses that corresponds to a 
particular value of prestrain. 

To ensure crack nucleation and growth at the nanometer scale, the prestress was varied 
between 3.5 and 4.25 GPa. If the system were uniaxially strained, these very high prestress 
values would have caused strong plasticity effects not related to the crack, such as spontaneous 
dislocation nucleation from the GBs (Yamakov et ah, 2002). Applying triaxial hydrostatic stress 
eliminates these undesirable plasticity effects. 

The crack in the system is nucleated by screening (preventing) the atomic interactions 
between atoms at both sides of the GB plane between Crystal I and Crystal II along a region of 
length l Q . As the crack grows and the crack opening becomes large enough to prevent interaction 
of atoms on the adjacent crack faces along the screened region, the previously screened atomic 
interactions are restored and the crack continues to evolve on its own. The crack starts growing if 
/ 0 is larger than the critical Griffith length L g defined when the energy spent to create the upper 
and lower crack surface 2y s minus the energy gained by destroying the GB y G B is equal to the 
released strain energy -dU/dl, per length /, 

2y s - y G B = -dU/dl. (2) 

An estimate of L g is made by calculating dU/dl as a function of a and /. This calculation is 
performed by using an anisotropic elastic finite element model of the elastic equivalent of the 
MD system, as will be discussed later. For the values of prestress applied in this study, a = 3.5, 
3.75, 4.0, and 4.25 GPa, the obtained Griffith lengths are L g = 6.08, 5.32, 4.71, and 4.21 nm, 
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respectively. The initial crack length is set equal to the crystallographic period in the x-direction, 
/ 0 = <7 7 1 0>a o = 5.7 nm, to reflect the existing periodicity of the equilibrated GB structure 
(Dahmen, 1990) induced by the lattice of the two joined crystals. This value of l 0 is larger than L g 
for most of the prestresses, except for o = 3.5 GPa, ensuring the initiation of crack growth for 
these prestresses. Crack growth also occurs for a = 3.5 GPa ( L g = 6.08 mn), because the use of a 
many-body potential causes all the atoms within the interaction range of the screened atoms to be 
also affected by the screening, effectively increasing the initial crack length. 

The identification of various structural defects including dislocations, twins, stacking faults, 
etc. appearing around the growing crack is important for understanding the mechanism of 
deformation that influence the CZM constitutive law. A procedure for atom identification based 
on the atom’s coordination number and on the common-neighbor-analysis (CNA) technique 
(Honeycutt and Andersen, 1987; Clarke and Jonsson, 1993) is used. The technique makes it 
possible to identify atoms in fee and hep states. Layers of hep atoms in an fee lattice are formed 
at stacking faults and twin boundaries (Weertman and Weertman, 1992) and can be used very 
successfully for visualizing the ongoing dislocation processes in fee crystals (Schiotz, 1999). 
Atoms that are not identified in an fee or hep state are considered to be in a non-crystalline state 
and indicate the presence of GBs or dislocation cores. In addition, atoms that have lost more than 
1/3 of their neighbors inside the interaction range of the potential are considered surface atoms. 
On average, 1 nm 2 of a flat {5 5 7} surface contains 16 surface atoms. The crack free surface S is 
estimated by counting the number of surface atoms, and when divided by 2 h (see Fig. 1, the 
prefactor 2 accounts for the two crack surfaces), conveniently gives an effective crack length /. 
Under this convention, the thermalized structure of a (l 10)299 symmetric tilt GB, shown in Fig. 
2, appears as quasi-periodic, with a regular pattern of hep atoms immersed in a disordered layer 
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with the presence of a few distributed vacancies, identified by the appearance of a few isolated 
surface atoms inside the GB region. As will be shown later in this paper, classifying atoms in this 
way presents a unique possibility to distinguish and quantify the various atomic processes 
occurring at the crack tip, and to identify their contribution to the decohesion law of the CZM 
elements. 

2.3 The finite element model 

A FE configuration, with dimensions scaled to reproduce the proportions of the MD system, 
is presented in Fig. 3. The FE calculations are anisotropic linear-elastic, carried out under 
displacement control in the x- and y- directions, and use generalized plane strain to reproduce the 
hydrostatic triaxial stress in the MD simulation. The use of periodic boundary conditions in the 
constant strain MD simulation, described previously, constrains the problem in the same manner 
as the generalized plane strain condition in continuum mechanics. The commercial software 
package ABAQUS (ABAQUS, 2004) is used. The model contains a combination of six-node 
triangular elements and eight-node quadrilateral elements that support generalized plane strain. 
The anisotropic elastic constants are obtained from the MD interatomic potential and 
transformed according to the crystallographic orientations of the atomic microstructure (see 
Appendix A). Higher order terms in the elastic constants are not included as they are not 
available from the published data for this potential (Mishin et ah, 1999) and are not easy to 
calculate. This limits the FE simulation to be linear elastic. 

The FE model contains a built-in lenticular slit, which simulates the MD approximation to a 
crack of varying relative length, 0.05 < l/L < 0.91, and is used to study the evolution of the 
system at different stages of the crack growth. To avoid stress singularities, the edges of the slit 
have a finite, but very small, initial radius X = 0.008/., which, when scaled to the MD 
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dimensions, corresponds to an initial crack opening of approximately 1 mn, i.e., slightly larger 
than the range of the interatomic MD potential. 

Two significant differences between the FE and MD models have to be outlined here. First, 
the FE calculations are static and linear-elastic. No dynamic or nonlinear effects are 
incorporated, and the results are relevant for static “cracks” under small strain when the 
nonlinear effects can be neglected. Second, there is no elastic equivalent of the GBs in the FE 
system. It is assumed that the difference between the elastic response of the GBs and the interior 
of the grains does not significantly alter the crack behavior because of the relatively small 
volume ratio of the GB in the system. Thus, the FE simulations used here play mainly a guiding 
role in setting up the parameters of the MD model and in understanding the MD results by 
distinguishing between plastic and elastic processes in the MD simulation. 

The first implementation of the FE model is to perform energetic analysis of the system in 
the absence of nonlinear elastic or plastic processes and to extract the Griffith length L g for the 
different loads of the MD system using Eq. (2). The elastic strain energy per unit thickness, E/h, 
of the FE system versus increasing crack length, 0 < / < 110 mn, is presented in Fig. 4. The 
values for / = 0 relates to a system with no crack, and is given for reference. Because the FE 
simulations are linear elastic, the strain energy can be easily scaled to the absolute dimensions of 
the MD system and to the four different applied prestresses a. Knowing y s and Ygb, extracted 
from the MD model and given in Section 2.2, one can plot their effective contributions to the 
crack energy as two straight lines: for the surface energy, £ slir f = 2y s /; and for the change of the 
GB energy, AEgb = -Ygb/, as shown in Fig. 4. Summing the elastic energy E e i ast (o,l), £ sur f, and 
AEqb gives the total energy balance of the system for each prestress condition, 

E t0 t(o,P) = Eeiast(a,/) + (2y s - Ygb)/ . (3) 
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For each a, E tat (o,l) has a maximum at the Griffith length / = L g {o) as shown in Fig. 4. The 
precise values of these maxima are obtained by using second order polynomial interpolation of 
the FE data points for E e \ ast (o,l) at / = 0, 6, and 10 mn in Fig. 4 for each prestress a. These results 
set the initial length, I 0 =5.7 mn, of the nucleated crack in the MD model, as discussed in the 
previous subsection. More importantly, the FE simulation shows the region for steady-state crack 
propagation. This steady-state region is where the decrease of Got(cx,/) with increasing / becomes 
almost constant. The steady-state region is revealed by drawing a straight, or steady-state, line ( s - 
s ) through EeiastfcT,/) for o= 4.25 GPa in Fig. 4. The line overlaps with Got in the interval 30 < / < 
80 mn indicating the steady-state region in Fig. 4. For L g < l < 30 mn, the crack propagation is 
expected to increase ( nucleation region in Fig. 4) because the slope of Got increases, approaching 
the steady-state value. For / > 80 mn, the crack starts to be influenced by the finite length of the 
system {finite-size effects region in Fig. 4), and Got increases above the steady-state value due to 
the increasing contribution of the crack’s periodic image. As the size of the steady-state region is 
a result of the system dimensions only, it remains the same for all the applied loads in the 
simulation. Thus, FE simulations help to estimate L g and to find the region of crack lengths 
where the crack is expected to grow under steady-state conditions, allowing for time independent 
statistics to be extracted for derivation of the CZM decohesion law. 

3. The mechanisms of crack propagation along the 299 grain boundary by MD 
simulation 

As explained in Section 2, the GB opens after 8 ps of screening of the atomic interactions in 
a region of 5.7 mn length along the middle of the GB between Crystal I and Crystal II (Fig. 1), 
and a crack starts to grow in both directions along the GB interface. Fig. 5 shows MD snapshots 
of cracks that have grown for four different initial hydrostatic prestresses: a = 3.5, 3.75, 4.0, and 
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4.25 GPa. In all cases, the crack growth is not symmetric in the +x and -x directions (as defined 
in Fig. 1) along the GB (® in Fig. 5). The growth in the -x direction has produced two twin 
patterns (© in Fig. 5) by emitting a series of twinning dislocations (© in Fig. 5) (Weertman and 
Weertman, 1992), and is almost symmetric about the grain boundary of the two adjacent crystals. 
The crack propagation in this direction is greatly reduced, when compared to the propagation in 
the +x direction, with energy expended through ductile blunting at the crack tip. The growth in 
the +x direction proceeds through a continuous process of void formation (© in Fig. 5). Because 
the crack tip emits very few dislocations (© and © in Fig. 5), crack propagation in the +x 
direction proceeds much faster than in the -x direction, and does so in an almost brittle fashion. 
After being emitted, some of the dislocations glide away from the crack tip and are absorbed by 
the next GB layer to form GB dislocations (© in Fig. 5). Secondary slip (© in Fig. 5) on a slip 
plane along the [l 1 2] direction (see Fig. 2) may also accommodate deformation non-parallel 

to the primary slip plane along the [l 1 2] direction, where twinning ® and primary slip 

dislocations © propagate. 

A detailed analysis of all of the atomistic processes, accompanying the crack growth and 
governing the decohesion of the GB interface follows. 

3.1 Twinning at the crack tip - crack propagation in the -x direction 

Fig. 6 contains a series of snapshots of the crack shown in Fig. 5(c) monitoring the formation 
of two symmetrical twins at the crack tip propagating in the negative, -x, direction in the system 
(as defined in Fig. 1). Although twinning is not common in Al, twin formation near a crack tip 
has been reported on several special occasions, both in MD simulations (Farkas et al., 2001; Hai 
and Tadmor, 2003) and in experiments (Pond and Garcia, 1982; Chen et al., 1999). Twinning has 
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been found to occur in special cases, when the crack propagation direction is suitably oriented 
with respect to the available slip planes in the crystal. In these references, the crack propagation 
is transgranular and the twins are formed along the most energetically favorable slip plane 
inclined with respect to the crack propagation direction. 

As shown in Fig. 6(a), the twinning process starts with a spontaneous emission of 1/6 

[1 1 2] Shockley partial or twinning dislocations © (see, also, © in Figs. 5(a-c)) (Weertman 

and Weertman, 1992) on the (1 1 1) planes in both crystals (see Fig. 2 for the crystallographic 
orientations in the model). Each dislocation creates a stacking fault ©, which is broadened into a 
nanotwin © by the subsequent emission of more twinning dislocations © of the same type (Figs. 
6(b-f)). Each new dislocation results in the propagation of the crack tip by one atomic layer. If 
the twinning process starts simultaneously on both, the -y face and the +y face at the -x crack 
tip, as shown in Fig. 6(a), the two twins grow symmetrically in Crystal I and Crystal II (Figs. 
6(b-f)). If the twinning processes do not start simultaneously due to some small perturbations in 
the simulation, as shown in Fig. 7(a) (to be discussed later), where the twinning in the -y 
direction © is initiated first, the two twins in Crystal I and Crystal II may be offset (Fig. 7(b)). 
As the crack grows, the symmetry is nearly restored, as can be seen in Fig. 5(b). Thus, the 
symmetry of the crystal lattice on both sides of the GB leads to a corresponding symmetry in the 
twin formation, which follows the symmetric orientations of the slip planes relative to the GB. 

Secondary slip inside the twins occurs in the fonn of partial dislocation emission identified 
by trailing stacking faults (© in Fig. 5(c) and Figs. 6(e,f)). This happens because the growing 
twin cannot accommodate all of the shear produced at the crack tip. The twinning changes the 

orientation of the (1 1 1) slip plane from ai ~ 9° to a new orientation a.i ~ 30° from the crack 
plane (Fig. 6(f)). The latter orientation is more favorable for dislocation emission because it 
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suggests a higher resolved shear stress. Thus, secondary slip is initiated on that plane, giving an 
additional degree of freedom, which makes accommodation of any deformation in the (x,y) plane 
near the crack tip possible. 

The active twinning through the spontaneous emission of twinning dislocations at the crack 
tip, together with the development of secondary slip inside the twins, makes the propagation in 
the -x direction of a ductile type and involves plastic mechanisms that expend considerable 
energy. The MD simulation results from this crack tip will later be used to construct the cohesive 
law for a ductile type CZM element. 

3.2 GB decohesion by cleavage - crack propagation in the +x direction 

The propagation by cleavage in the +x direction for the crack under a 3.75 GPa prestress is 
monitored in Fig. 7. First, the disorder of the GB region in front of the crack tip is increased (Fig. 
7(a)) until a void appears (© in Fig. 7(b)). The void grows (© in Fig. 7(c)) and eventually joins 
the main crack, increasing its length (Fig. 7(d)). Then, another void forms (© in Fig. 7(d)) and 
the process repeats, becoming the primary mechanism for crack propagation. Figure 7(e) shows 
the formation of two successive voids in front of the crack tip that are about to join to form an 
increment of crack growth (Fig. 7(f)). A similar mechanism of fracture in an alloy system has 
been reported in MD simulations of intergranular crack propagation along a high-angle GB in 
NiAl (Farkas, 2000a; 2000b). 

Interestingly, dislocations are not emitted until the final stage of growth (Fig. 7(f)), when the 
crack is blunted by a spontaneous emission of dislocations (© and © in Fig. 7(f)). The 
spontaneous dislocation emission is believed to be caused by a dynamic instability at the crack 
tip that is known to happen at a propagation speed of approximately 1/3 the Rayleigh speed of 
sound cr, and leads to a dynamic brittle-to-ductile transition in the crack propagation (Abraham 
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et. ah, 1994; Abraham, 2001; Gao, 1996; Yamakov et ah, 2005). The appearance of the dynamic 
instability in the theoretical model of Barber et al. (1989), used as a base for the system 
simulated here (discussed in Sec. 2.1), was predicted by Ching (1994). There, it was stated that 
steady-state crack propagation above the critical speed is unstable. The MD simulation, in the 
present work, supports this statement. Additional details are given by Yamakov et al. (2005), 
wherein the dynamic aspects of crack propagation are discussed. 

The lack of dislocation emission from the crack tip propagating in the +x direction localizes 
the damage zone along the GB interface, where the microvoids are nucleated (Fig. 7(b-e)). The 
+x propagation is much more brittle and propagates much further along the GB than is observed 
in the —x direction. The results from this crack tip will later be used to construct the cohesive law 
for a brittle type CZM element. 

3.3 Rice and Tadmor-Hai criteria for twinning vs. cleavage at the crack tip 

The ability of GBs in fee metals to produce asymmetric crack growth was first found in an 
MD simulation of intergranular fracture in Cu by Cleri et al. (1999). In their work, a crack 
propagated in the interface plane of the symmetric tilt (22 1)/(22 1) grain boundary. The crack 

advanced by brittle fracture along the [1 1 4] direction and was blunted by dislocation 

emission along the opposite [1 1 4] direction. The difference in behavior at the two crack tips 
was attributed to the orientation of the slip planes relative to the GB. A slip plane inclined at 
angle 0 to the GB makes the angle 0 to the propagation direction of one crack tip, and the angle 
Tt-0 with the propagation direction of the opposing tip. Thus, for certain GBs, the Rice criterion 
(Rice, 1992) for dislocation nucleation versus cleavage might be satisfied for the crack tip at 
angle 0, but not for the crack tip at an angle rc-0, leading to asymmetric crack propagation. The 
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crack grows primarily in the direction of the jt-0 tip since the brittle propagation takes much less 
energy. Note that, for transgranular fracture, every direction in the fee lattice has inverse 
symmetry, so such asymmetry in crack propagation is not expected to appear. 

In the crystallographic model, presented in Fig. 2, the propagation direction [7 7 10] 
along the GB plane makes an angle 0=79.98° (see Fig. 6(f)) with the twinning direction 

[1 1 2] in both crystals, giving cos0 = 1 / a/ 33 . The Rice criterion for dislocation emission for 
the two propagation directions with +cos0 (-x direction) and cos(jt-0) = -cos0 (+x direction) 
can be written as 

Y,„ < (l — cos0)sin~ 0 

2 Y. S " Ycb 8 

where the dislocations at the crack tip are edge twinning dislocations (shown as (D in Fig. 6) with 
their Burgers vector coinciding with the twinning direction and lying in the (110) plane of the 
model. Using the values for y us , y s , and yoB given in Sec. 2.2, the left side of Eq. (4) equals 
0.127±0.01. The right side is determined to be 0.142 and 0.100 for the -x and +x propagation 
directions, respectively. Thus, the Rice criterion predicts that the crack tip should propagate in a 
ductile manner in the -x direction, and in a brittle manner in the +x direction, in agreement with 
the behavior seen in Fig. 5. 

While the Rice criterion explains the ductile-brittle asymmetry of the crack found in the 
simulations (Figs. 6 and 7), an additional consideration is needed to explain the preference for 
twinning rather than the emission of slip dislocations in the -x direction. Tadmor and Hai (2003) 
have recently derived a Pierels based criterion for the onset of deformation twinning at a crack 
tip in fee metals by comparing the energies for nucleating a twinning or slip dislocation. This 
criterion introduces a so-called unstable twinning energy y ut , which, by analogy with the unstable 
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stacking-fault energy y us , gives the energetic barrier that must be exceeded to form a twin from 
an existing stacking-fault. The criterion is given by the parameter T as 
f<l slip dislocation 

(5) 

[> 1 deformation twin 



where X is a factor that is dependent on, both, the y s f / y us ratio, and the orientation of the slip 
system of the twinning dislocations with respect to the crack plane (Tadmor and Hai, 2003). For 
the twinning tip in this simulation X=1.51 (as applied from Tadmor and Hai, 2003), and 

= 0.9, which results in T =1.36. This suggests that twinning, rather than slip dislocation 



emission, will be the preferred deformation mechanism. 

The brittle-ductile analysis based on the Rice criterion is for static or sufficiently slow crack 
growth. The criterion does not account for dynamic effects occurring at the crack tip propagating 
in the +x direction (Yamakov et ah, 2005) and cannot explain the observed dislocation burst in 
Fig. 7(f) (discussed in Section 3.2). Crack propagation is expected to be brittle throughout. 
Nevertheless, apart from these dynamic effects, there is good agreement between the simulation 
results and the two conditions expressed through Eqs. (4) and (5). This allows the prediction of 
the type of decohesion by knowing the crystallographic orientations of the grains and a few 
material parameters, such as the surface, GB, stacking-fault, and twinning energies. 


4. Plastic contribution to the stress field near the crack: MD - FE comparison 

The plastic processes at the crack tips, including twinning and dislocation emission, have a 
pronounced effect on the stress distribution near the growing crack. This effect is best revealed 
by a comparison between the stress distributions obtained from the MD and the linear elastic FE 
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simulations for the models presented in Section 2. While the MD simulation considers the 
material structure at the atomic level and intrinsically incorporates all the plastic processes 
together with the elastic response of the bicrystal, the FE model assumes an elastically equivalent 
system and gives only the elastic solution. Another aspect of the FE simulation is that it gives a 
static solution when the system is in elastic equilibrium. By contrast, the MD simulation 
describes a continuously evolving system with all dynamic effects present. To avoid the 
implications of the dynamic effects in the comparison between FE and MD results, the 
corresponding MD stress distribution is calculated after the atomistic system reaches elasto- 
plastic equilibrium and the crack stops growing. The dynamic solution for the elastic stress as a 
function of the crack velocity in the steady-state analytical model (Barber et ah, 1989) is given 
by Ching (1994). Ching’s solution is not studied here because extracting the non-equilibrium 
dynamic stress from MD is not a trivial task and goes beyond the scope of this work. 

To obtain the continuum stress distribution from the MD simulation at elasto-plastic 
equilibrium, the virial formulation for the local stress is used, defined as (Cormier et ah, 2001) 


, 1 VI (0 (>) (>) 1 V ^£7 r a J), i' J) 

'«r^ " r «j> 

'EQ\ j ' 


where nr ' and v (!) are the mass and velocity of the i - th atom, respectively, and r u ' n = r 1 1 - r {n is 


,U,j) _ Jl) Jf) 


the vector between the positions of atoms (/) and (j) with a, (3=1, 2, 3 for the x, y, and z Cartesian 
vector components. The first sum is taken over atoms (/) in a volume Q over which the stress is 
calculated, and the second sum is taken over atoms (j), which are in the interaction range of the 


atom (/). The derivative of the interaction potential 


is a generalization of the pair 


interatomic force between atoms (/) and (j) in the case of a many body potential. The volume Q 


in this study is of size 6a 0 x 6 a 0 x h in the x-, y-, and z- directions, respectively, which defines 
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£2=16.89 nrn 3 containing 864 atoms. This volume is large enough to ensure a good convergence 
of the calculated stress, but small enough to allow for a sufficiently high resolution of the stress 
distribution map. Additional time averaging of the stress is performed over 16 ps of simulation 
time, which is sufficient to smooth out the noise from the phonon waves in the simulation. 

The virial theorem, which is the basis for calculating the local virial stress (Eq. (6)), provides 
the oldest and most frequently used expression for relating forces and motion within an atomic 
system to a continuum stress. The virial stress is defined through the local momentum flux 
carried by the particles in a small volume element (Lutsko, 1988), rather than as a force acting 
over a small surface element, i.e., the definition of the Cauchy stress used in continuum 
mechanics. Because the Cauchy stress assumes continuous structure of matter, the two 
definitions increasingly diverge at the atomic scale. However, in the limit of time and volume 
averages at equilibrium, the virial stress coincides with the Cauchy stress (Zimmerman et ah, 
2002 ). 

The comparison between the stress distributions obtained from the FE and MD simulations 
for two cases of cracks are given in Figs. 8 and 9 in the form of two-dimensional (x-v) stress- 
maps of the system created for each of the three stress components acting in the (x-y) plane: o xx , 
Oyy, and o xy . Fig. 8 presents the case of the small prestress of 3.5 GPa, which produces a 10 nm 
crack at equilibrium in the MD simulation (seen in Fig. 5(a)). The size of the crack is sufficiently 
small compared to the system size, and the elastic stress field calculated by the FE simulation 
(Figs. 8(a-c)) does not experience edge effects from the boundary conditions and thus resembles 
an infinite plate. The stress distribution for the brittle crack tip (the circled area in Figs. 8(d-f)) 
obtained by the MD simulation shows very good quantitative similarity to the FE result. This 
similarity suggests that the crack propagation in the MD system in the +x direction is indeed 
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brittle with essentially no plasticity. However, for the twinning tip there is a significant 
difference between the MD and FE results. The FE stress distribution is symmetric for the two 
crack tips, as they are elastically equivalent in the FE model (Figs. 8(a-c)). In contrast, twinning 
occurs in the MD system at the tip propagating in the -x direction and significantly alters the 
stress field. The o xx and o xy stress components are largely relieved, while the o yy stress is 
distributed along the twin boundaries (© in Fig. 8(e)) and relieved at the crack tip. This 
reduction in local stress state explains the very slow crack propagation in the -x direction of the 
MD model. 

The GB layer, present in the MD model, but not in the FE model, increases the o xx stress (® 
in Figs. 8(d) and 9(d) corresponding to ® in Figs. 5(a,b)), and suggests a stiffer GB layer. This 
stiffening effect should be taken with caution as the virial stress calculation in the MD simulation 
(Eq. (6)), as discussed previously, may give erroneous results when computed over small 
disordered atomic domains containing structural defects (e.g., GBs and surfaces, dislocation 
cores, etc., Zimmerman et ah, 2002). It has been suggested that the highly constrained 
thermodynamic state of GBs in multilayer systems may appear stiffer than the crystal lattice in 
certain directions (Wolf and Jaszczak, 1993). 

At a larger prestress of 3.75 GPa the crack in the MD model reached almost 40 nm in length. 
The stress fields, presented in Fig. 9 (corresponding to the MD snapshot in Fig. 5(b)), are 
strongly affected by the system boundaries at y = ±W/2, as compared with the stress fields in Fig. 
8. In the MD case (see Figs. 9(d-f) related to the snapshot in Fig. 5(b)), the stress field away from 
the crack tips is similar to the corresponding FE stress field, while the near-tip stress distribution 
is strongly affected by the plastic mechanisms at the crack tips. The contribution to the stress 
from the developed twinning on the -x side of the crack (© in Figs. 5(b) and 9(e)) and the few 
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dislocations remaining in equilibrium with the crack tip on the +x side (© in Figs. 9(d-f)) can be 
easily identified. Both processes effectively relieve the stress around the crack tips as compared 
to the FE elastic solution (Figs. 9(a-c)) and lead the MD system into elasto-plastic equilibrium. 

The stress-field maps discussed in this section are the basis for continuum interpretation of 
the atomistic MD simulations. They show that the stress around a crack tip, computed in the MD 
model, can be related to a continuum stress field. This makes it possible to use the MD 
simulations of crack propagation to extract the cohesive zone law for a grain-scale continuum 
model. 

5. Defining a Traction-Displacement Relationship from Molecular-Dynamics 

Grain-scale simulations that use cohesive zone models to study fracture (such as those 
presented by Zavattieri et al. (2001; 2003), Iesulauro et al. (2002), and Wei and Anand (2004)) 
typically use heuristically derived relationships and input values to define the CZM. In such an 
approach, the values at the microstructural scale are based on macroscopic parameters such as 
fracture toughness. Thus, there is no physical substantiation for the quantities that are input to the 
analysis. Conversely, a cohesive formulation can be part of an effective physics-based approach 
if constitutive parameters are used that exploit the similitude between atomistic simulation and 
continuum finite element results. The model can be developed from results of the MD analyses 
allowing the physical insight of MD to be embedded in the more computationally efficient FE 
models. 

The traction-displacement function of a CZM describes how the traction x, developed in the 
process zone near the crack, depends on the crack opening X. In general, the traction- 
displacement function x(X) considers both normal and tangential crack surface components of the 
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traction. In the case of a normal opening (mode I) crack under hydrostatic load, only the normal 
component of the traction is important. In the MD simulation, the tensile a s vv component of the 
stress near the debonding GB interface can be used as an equivalent of the normal traction x n . 
The superscript V indicates the stress near the GB interface. 

To obtain a s vv from MD, the simulation box is divided into 2N[ + \ slices of thickness 6 (Fig. 
10), so that (2/V/ + 1 )5 = L. At the beginning of the simulation, before the crack is nucleated, a 
three-dimensional rectangular volume element (called a Cohesive-Zone- Volume-Element, 
CZVE) is defined at the place where each slice crosses the GB plane at y = 0 (illustrated by the 
dotted regions in Fig. 10). The CZVE has a length 8* along x and extends to y = ±5 V on both sides 
of the GB interface. Its thickness equals the system thickness h. When the propagating crack 
passes through a CZVE during the simulation and the CZVE starts to open, its shape deforms 
(see the dotted elements at the crack tips in Fig. 10). To be able to keep track of the CZVE 
volume, the atoms that have belonged initially to the CZVE are marked, i.e., their identification 
numbers are stored for each CZVE, and are further used to identify the CZVE. Eq. (6), with the 
“/’’-index going over the marked atoms of a CZVE, and Q equal to the volume of the CZVE, is 
used to calculate c f vv (x p ,t) as a function of the position x p of the p - th CZVE along the GB 
interface at time t during the simulation. Here, p G [0, Nl] for the crack tip propagating in the +x 
direction, and p G [0, -Nl] for the crack tip propagating in the -x direction. At the same instant t, 
the crack opening X(x p ,/) at the same position x p , is also estimated. Thus, every estimate of 
c f n (x p ,t) corresponds to one estimate for k(x p , /), allowing the construction of the function 
c f V y(k(x p ,t)) for each CZVE. 

An example of estimating c i s vy {x p ,t) and ’k(x p ,t) for the case of a crack growing for t = 123 ps 
in a system prestressed at 4.25 GPa is given in Fig. 11. The stress and opening profiles (o' s >> (x) 
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and X(x) curves, respectively) are taken from a set of CZVEs placed in a line along the 
debonding GB. A fit to the continuum expression for the stress at a crack tip 


o' W = 


2k(x - x ) 


is imposed on the stress profile ahead of the brittle tip, propagating in the +x direction (see Fig. 
11). The fit shows a 1 Nr type dependence of the stress profile and Ki is detennined to be 0.36 
MPa.m . Such dependence is not recovered for the ductile tip, propagating in the -x direction in 
Fig. 1 1, because the stress is continuously relieved by the twinning process. 

The c fyy(x) and X(x) profiles in Fig. 11 are for a single instant of time of the crack 
propagation. If the whole simulation of the interface debonding is divided into N, equal intervals 
of time t q (gE[0, A,]), many such profiles can be taken of many CZVEs placed along the GB. 
When plotted in a vs. X plot, each (c r\ y (x p ,t q ), X(x /; , t q )) couple represents a point c f vy (k(x p ,t q )). 
After sorting these data points in an ascending order on X, so that X/ < X/+i, and taking a moving 
average (or a consecutive mean): 


, M 


in which the results are averaged over M points backward and M points forward from X/, a 
construction of a statistically representative traction-displacement function x(X) can be made. 

The opening of the crack X(x,t) as a function of x can be estimated in various ways. A 
convenient method is to use the previously described slicing of the simulation system along the 
x-axis (Fig. 10). The separation between atoms on either side of the GB interface increases as the 
crack passes through a slice, creating an opening of size X in the middle of the slice (the dashed 
slice in Fig. 10 is given as an example). The opening changes the slice’s gyration radius R g (or 
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the angular momentum), and one can estimate the size of the opening by calculating R g . 
Assuming a uniform mass distribution, the continuum expression of R g as a function of X is 


W 12 

}y 2 dy 

= — (w : + wi + V). 

5“y 12 

X/2 


( 9 ) 


Alternatively, expressed through the atomic coordinates, the discrete formula for R g is: 


R 


2 

g 


N 


2 ( yi-ycf 




( 10 ) 


where y,- is the y-coordinate of atom i, and the sum is over all N atoms in the slice. y c = 



is the y-coordinate of the mass center of the slice. Calculating R g from Eq. (10) and substituting it 
in Eq. (9), the crack opening at the slice is the positive root of X from Eq. (9). 

Though the introduced CZVE is used to extract the behavior for the CZM element in the FE 
model, it should not be considered an MD equivalent of a CZM element. While the CZM 
elements are strictly surface elements in a 3D finite element model, the CZVE encompasses a 
volume at the surfaces on either side of the crack, large enough to allow for the convergence in 
the stress estimate in (Eq. (6). Each CZVE along the GB interface in the MD simulation has its 
individual behavior, depending on the local structural irregularities at the atomic level, while in a 
FE simulation the CZM elements for one type of interface are all identical. The relation between 
the CZVE and the CZM is made through statistical averaging (Eq. (8)). The assumption is made 
that under the steady-state conditions of crack propagation, the statistical average over many 
CZVEs scanned at different instants of time will produce a statistically unique traction- 
displacement curve. Assuming a Gaussian distribution of c?yy(k(x,t)) from each individual 
measurement of each CZVE, the dispersion around the average traction-displacement curve will 
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decrease proportionally to the square root increase of the number of atoms belonging to the 
CZVE. Thus, a larger volume of the CZVE would better average out the statistical fluctuations. 
In MD simulations, this volume is limited by the small size of the damage zone around the crack 
tip and by the large stress gradients inside it. This makes the dispersion large, and a sufficient 
number of measurements is required to produce a reliable traction-displacement function. 

The size of the CZVE is defined by 6. v = b y = b = 1/3 [7 7 10]a o ~ 1.9 mn (see Fig. 10) for this 
simulation, which gives a volume of 2hb 2 = 20.7 mn 3 containing 1245 atoms that contribute to 
the stress-displacement response of one CZVE. For the simulated system in Fig. 1, 2 Nl + 1 = 63 
CZVEs fit along the 120 mn GB interface. During the simulation, cf w (X) state for each CZVE is 
scanned every 4 ps for a period of 200 ps crack propagation time. This gives 3150 points for 
c fyy(X) dependence from each simulation run. These points provide adequate statistics to extract 
the traction-displacement curve through suitable averaging. Because the two crack tips propagate 
in a very different manner, as discussed in Section 3, the statistical averaging (Eq. (8)) must be 
performed separately for these two directions. 

The values of cfy/Z) for the brittle tip prestressed at 4.25 GPa are given in Fig. 12(a). At 52 
mn of propagation (Fig. 5(d)), 25 CZVEs have experienced a complete transition from a fully 
closed to a fully opened state. The results from 1575 calculations (along 0 < x < C/2), presented 
as points in Fig. 12(a), gives sufficient statistics to extract a reliable cr s yy (X) curve using Eq. (8). 
The averaging interval is set to M= 100, which is empirically found to give a sufficiently smooth 
curve. The curve reproduces closely a bilinear type of constitutive relation for the CZM element 
and is similar to the one suggested by Camacho and Oritz (1996) for brittle materials. The 
bilinear types of CZM models were used in several recent FE simulations of brittle fracture 
during multi-axial dynamic loading of ceramic microstructures (Zavattieri et al., 2001; 2003). 
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The traction-displacement curves for 3.75 and 4.00 GPa, obtained in the same way, also given in 
Fig. 12(a) (the data points, over which the traction-displacement curves are averaged are omitted 
for clarity), show that varying the prestress of the system in this range does not change the 
essential character of the averaged o' s >y (X) dependence. The peak stress o p shows a slight 
systematic decrease with decreasing load, and is probably an effect of the propagation speed 
decrease but remains approximately o p ~ 5.0 GPa at k p ~ 0.4 mn. It should be noted that the 
value for is taken relative to the initial state at an applied hydrostatic loading of 3.75 to 4.25 
GPa. In reality, o' s >y (X) for K < k n should follow the elastic response of the GB before debonding. 
The full opening, or full debonding, of the interface happens at k 0 ~ 2.5 mn, when & yy becomes 
zero. 

The corresponding data o' s >> .(X) for the ductile tip, propagating in the —x direction where 
twinning occurs, is presented in Fig. 12(b). Again, only the data points for 4.25 GPa prestress are 
given for clarity, while the averaged traction-displacement curves are shown for the three loads 
of 3.75, 4.0 and 4.25 GPa, as in Fig. 12(a). In contrast with the behavior in the +x direction, a 
stronger influence of the prestress on the traction-displacement curves is observed for the —x 
direction. Most pronounced is this influence for k 0 , which changes from h ( { 1 1 = 2.7 mn to A (l <3) = 
4.3 mn as the prestress is increased from 3.75 to 4.25 GPa. The peak stress o p shows again the 
slight systematic decrease with decreasing load, which was observed in Fig. 12(a), accompanied 
by a shift in X p . 

The comparison between the extracted traction-displacement relationships for the ductile and 
the brittle crack tips is presented in Fig. 13. The case of 4.25 GPa prestress is chosen because 
there, the crack has achieved its largest growth in both directions (Fig. 5(d)) and has produced 
the most data points for CZVE statistics. Figure 13 reveals how the different atomistic processes, 
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taking place at the two crack tips and discussed in detail in Sec. 3, are affecting the decohesion 
law of the interface debonding. This comparison is also very informative on the significance of 
the form of the traction-displacement curve and how this form reflects the underlying physical 
processes as revealed by the presented atomistic MD simulations. The ductile traction- 
displacement curve is substantially more extended towards larger debonding distances X 0 , 
compared to the case of almost brittle debonding (k 0 ' /f > 'k 0 br , see Fig. 13). The corresponding 
peak stress o p is also lower for the ductile case, which is possibly a result of the emission of 
twinning dislocations relieving the stress at the crack tip. The lower o p and the larger A 0 during 
twinning produces a more extended debonding region for the ductile traction-displacement 
curve, compared to the steeper curve in the case of brittle fracture. This difference justifies the 
use of a trapezoidal type of curve with a plastic straining region and reduced peak stress for the 
case of elasto-plastic fracture (Tvergaard and Hutchinson, 1996) and a bilinear curve for brittle 
fracture (Camacho and Ortiz, 1996). Despite a lower peak stress, the traction-displacement curve 
with twinning has about 50% larger area, meaning a larger work of separation T defined as 
(Tvergaard and Hutchinson, 1992) 

di) 

The larger work of separation also explains the smaller propagation distance in the ductile (— jc) 
direction, compared to the brittle (+x) direction in the MD simulation (see Fig. 5(d) and compare 
also Figs. 6(a-f) with Figs. 7(a-f) discussed in Section 3). 

6. Conclusions 

A methodology is detailed for extracting the decohesion law for interface debonding by 
introducing a Cohesive-Zone- Volume-Element using molecular-dynamics simulations. The 
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methodology is applied for the process of debonding of a high-energy 299 GB in A1 through 
crack propagation under hydrostatic loading conditions. The work reveals the atomic 
mechanisms in the damage zone near the crack tips during intergranular fracture. Here, the crack 
propagation has been shown to proceed in a different manner in the two opposite directions 
along the GB interface. While in one direction the crack progresses in a brittle manner, the 
propagation in the opposite direction is of a ductile type, characterized by the presence of plastic 
mechanisms including twinning at the crack tip. The difference in the mechanism is due to the 
inclination of the slip planes to the GB interface making different angles with respect to the two 
opposite propagation directions. 

Applying a statistical procedure, the decohesion laws in terms of traction-displacement 
curves for the brittle and ductile crack propagations are extracted from the atomic forces and 
motions around the crack tips. The comparison between these curves reveals the influence of the 
plastic mechanisms on the behavior of the CZM element. The appearance of plastic mechanisms 
extends the traction-displacement curve towards larger openings, while reducing the peak 
opening traction. This results in an increase of the area below the curve, as the work of 
decohesion increases in the presence of plastic processes. This justifies the use of an empirically 
derived trapezoidal type of decohesion law for the case of solids with enhanced plasticity and a 
bilinear law for britle materials. 

The brittle vs. ductile propagation simulated in both directions along the GB interface is 
found to agree well with the predicting criteria of Rice and Tadmor-Hai. This good compliance 
with the two criteria allows the selection of the proper cohesive zone model in a finite element 
simulation by knowing the crystallographic orientations of the grains in the micro structure and a 
few material parameters. The parameters, including surface, GB, stacking-fault, and twinning 
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energy, are readily available experimentally. Others, such as the unstable stacking-fault energy 
and the unstable twinning energy, need more precise estimates and may require the use of first- 
principal calculations. 
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Appendix A 

Molecular-dynamics (MD) monocrystal simulations were used to estimate the stiffness 
coefficients of an aluminum single crystal at 100 K in the basic orientation: (x:[l 0 0]; v:[0 1 0]; 
z :[ 0 0 1]). The values obtained are (in GPa): Cn = 109.1; C12 = 58.5; C 44 = 33.0. While C11 and 
C12 are slightly lower than the zero temperature value for the potential used (Mishin et ah, 1999), 
C 44 is slightly higher. This increases slightly the anisotropy parameter Ca = 2Cwl(C\\-C\?) = 1.30 
compared to the zero temperature Ca (T = 0 K) = 1.21. The transformation of the stiffness 
coefficients to the ^99 crystal orientations, given in Section 2.2, is performed according to the 
shortened for cubic crystals procedure (Wortman and Evans, 1965), which gives the following 
stiffness matrix 

118.7 52.7 54.7 0 0 ±2.0 

118.8 54.6 0 0 ±1.9 

116.8 0 0 +3.8 

29.1 +3.8 0 

sym. 29.2 0 

27.2 

where the prime indicates transformed elastic constants with C',y = C'p. The “±” sign at some of 
the off-diagonal values indicate change of sign when going from Crystal I orientation to Crystal 
II orientation. Due to the small anisotropy, the stiffness matrix in Eq. (A.l) almost preserves its 
cubic symmetry with C'n ~ C' 22 ~ C' 33 ~ 1 18, C'i 2 ~ C'i 3 ~ C' 23 ~ 54, C' 44 ~ C' 55 ~ C' 66 ~ 29, and 
the remaining C',y ~ 0. Using these approximate values in the FE model makes the two 
crystallographic orientations on both sides of the GBs elastically equivalent. This significantly 
simplifies the problem from a slit in a multilayer system to a slit in a single cubic crystal with the 
GBs neglected (see Fig. 3). 


(A.l) 
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Figures: 



Fig. 1 . Initial atomistic snapshot of the molecular-dynamics (MD) system before crack initiation, 
representing the simulation set-up, explained in the text. Size dimensions are in mn. Grain 
boundaries (GB) are shown as parallel lines of dark atoms separating the crystalline phases of 
Crystal I, Crystal II, and Absorbing Layers I and II as indicated. 
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Fig. 2. Atomistic snapshot giving the crystallography and structure of the GB interface. Common 
neighbor analysis (CNA) (Honeycutt and Andersen, 1987; Clarke and Jonsson, 1993) is used to 
identify atoms in different crystallographic states: fee (small dots), hep (triangles), and non- 
crystalline atoms (large dots). Atoms with more than 1/3 of their nearest neighbors missing are 
identified as surface atoms (squares), indicating existing vacancies in the GB. The length scale is 
in units of the lattice constant of Al, a 0 = 0.405 mn (Mishin et al., 1999). 
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Fig. 3. The mesh of the finite element model with a built in lenticular crack. The relative crack 
length, l/L, varied from 0.05 to 0.91, corresponding to the absolute crack length from 6 to 110 
mn in the MD system, shown in Fig. 1. Generalized plane strain in the z-direction is used to 
emulate the hydrostatic loading conditions in the MD system, shown in Fig. 1. 
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Fig. 4. The strain energy per unit thickness, E/h, of the FE system with increasing crack length /. 
The two straight, dashed, and dot-dashed lines, E sur f= 2y s / and A E G b = -y gbU show the effective 
increase of the surface energy and the decrease of the GB energy with increasing /. The position 
of the maximum of E tot = E e i ast + E sur f + &Egb indicates the Griffith length L g , as shown for the 
four applied prestresses of 4.25, 4.0, 3.75, and 3.5 GPa, starting from the top down. 
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Fig. 5. (In color) MD snapshots of cracks, which have propagated in the MD system, shown in 
Fig. 1, prestressed at four different initial hydrostatic stresses: a = 3.5 (a), 3.75 (b), 4.0 (c) and 
4.25 GPa (d). As in Fig. 2, CNA is used to identify atoms in different crystallographic states: fee 
(in gray), hep (in red), and non-crystalline atoms (in blue). Atoms with more than 1/3 of their 
nearest neighbors missing are identified as surface atoms (in green). Thus, a number of different 
formations are indicated in the figure as follows: © - GB interface; © - twin boundary; © - core 
of a partial or twinning dislocation; ® - nanovoid at the crack tip; © - slip dislocation; © - GB 
dislocation; and © - secondary slip. 
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Fig. 6. A series of snapshots monitoring the development of a symmetrical deformation twinning 
at the -x crack tip in the MD system prestressed at 4.0 GPa hydrostatic load. The snapshots are 
taken at various times from the crack initiation, t, given in ps. The length scale is in units of the 
lattice constant of Al, a 0 = 0.405 mn. Identified as in Fig. 5, the small dots indicate fee atoms; the 
non-crystalline atoms are shown in black, while hep and surface atoms are shown in grey. 
Formations ©, ©, ©, and © correspond to the types indicated in Fig. 5. 
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Fig. 7. A series of snapshots monitoring the propagation of the crack tip in the +x direction in the 
MD system prestressed at 3.75 GPa hydrostatic load. The snapshots are taken at various times 
from the crack initiation, t, given in ps. The length scale is in units of the lattice constant of Al, a 0 
= 0.405 mn. The atoms of different types are identified as in Fig. 6. Formations ©-© correspond 
to the types indicated in Fig. 5. 
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Fig. 8. (In color) Stress contours from the FEM model (a-c) (corresponding to Fig. 3) and the 
corresponding stress maps (d-f) from the MD model (corresponding to Fig. 5(a) at 3.5 GPa 
prestress) for o xx , o yy , and o xy stress components. In (a), (b), (d), and (e), the stress in blue is 
defined as the stress in tension. In (c) and (f), positive and negative shear corresponds to shear 
directions relative to the GB interface as shown at the two sides of the o xy stress indicator. 
Formations © and © correspond to those indicated in Fig. 5(a). 
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Fig. 9. (In color) Stress contours from the FEM model (a-c) and the corresponding 
(d-f) from the MD model (corresponding to Fig. 5(b) at 3.75 GPa prestress) for o xx , 
stress components. The stress-level colors are defined in the same way as in Fig. 8. 
©-© correspond to those indicated in Fig. 5(b). 
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Fig. 10. Schematic diagram of the slicing of the system volume in the MD simulation and 
defining the representative regions for extracting the parameters for the cohesive-zone interface 
elements in a continuum simulation. 
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Fig. 11. Stress and opening profiles extracted along the crack growing for 123 ps in a system 
prestressed at 4.25 GPa hydrostatic load. The corresponding snapshot of the crack is shown at 
the bottom. 
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Fig. 12. Surface stress vs. crack opening curves a s yy (k) characterizing the propagation of the 
cleavage tip in the +x direction (a) and in the -x direction (b) for three preloads. The dots are the 
individual measurements for 4.25 GPa hydrostatic load simulation, from which an average was 
taken to produce the corresponding a s yy (k) curve (see text). The curves for 3.75 and 4.0 GPa are 
produced in a similar way. 
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Fig. 13. Traction-displacement relationship a\ y (A) and the corresponding work of separation 
T(a) for the case of 4.25 GPa prestress compared for the brittle (+x) and the twinning (-x) tip. 
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